analytic Continuation of Quantum Monte Carlo Data: 
Optimal Stochastic Regularization Approach 



I. S. Krivenko and A. N. RubtscnQ 

Department of Physics, Moscow State University, 119992 Moscow, Russia 
(Dated: February 6, 2008) 

A new algorithm for analytic continuation of quantum Monte Carlo (QMC) data from the Mat- 
subara domain to real frequencies is proposed. Unlike the widely used maximum-entropy (MaxEnt) 
procedure, our method is linear with respect to the input data and can therefore be applied to 
off-diagonal components of Green's function, or to the self-energy function. The latter possibility 
is used to analyse QMC result for the half-filled Hubbard model on the Bethe lattice. Our method 
resolves the shoulders near the Hubbard band onsets, whereas the MaxEnt estimation does not. 
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I. INTRODUCTION 

Modern theoretical physics analysis of strongly cor- 
related systems deals with a wide range of numerical 
methods, because no serious progress has yet been made 
in formulating a regular analytic description. In par- 
ticular, progress has been achieved in Dynamical Mean 
Field Theory (DMFT)i^, as well as in its extensions 
and generalizations^. In the DMFT approach strongly- 
correlated systems are reduced to an effective impurity 
problem. Virtually, a single atom or small cluster is re- 
moved from the lattice but placed into a Gaussian bath 
defined in a self-consistent way. DMFT requires a so- 
called solver algorithm, which is intended to give an ap- 
proximate evaluation of Green's functions for an effective 
impurity model. For the most cases, the Quantum Monte 
Carlo (QMC) class of solvers is used^. These algorithms 
allow the calculation of electronic Green's functions for 
many models lying far beyond the regime where pertur- 
bation theory is valid. It is important to note that the 
output of QMC calculation belongs to the imaginary- 
time domain. This information is enough to evaluate 
DMFT loops, since DMFT equations are written in terms 
of fermionic Matsubara frequencies w n = (2n + 1)tv//3, 
where (3 is the inverse temperature. 

On the other hand, experimental data is, of course, ob- 
tained in the real-time or real-frequency domain. Conse- 
quently, it is evidently quite a difficult problem to extract 
physical information from a numerical data set, because 
analytic continuation is needed. 

In the canonical problem one has to obtain the spec- 
tral density function A(uj) from a noisy set of values for 
thermal Green's function G(t) or its Fourier transform 
Q{iuj n ) — J Q dr eT %UJnT Q(r), produced in a QMC simula- 
tion. The spectral density is proportional to imaginary 
part of retarded Green's function, from the definition: 

4(w) = --ImG» (1) 

7T 

Because of the analyticity of the retarded Green's func- 
tion in the upper complex plane of frequency the follow- 



ing is true: 
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Substituting uj = iu n one can derive the integral equa- 
tion 

= J du^pAV) (3) 

— oo 

which is a particular formulation of the analytic contin- 
uation problem. The kernel of this integral equation is 
exponentially small at large positive and negative fre- 
quencies, so tiny variations in Q(t) may correspond to a 
very strong change of the spectrum at those frequencies. 
Uncertainty in Q(r) is unavoidable due to the stochastic 
nature of QMC algorithms, so the problem of numerical 
analytic continuation is extremely ill-posed, similarly to 
a numerical inversion of a Laplace transform. 

The earliest attempts to solve the problem were based 
on Pade approximation method^. For a lot of cases, 
modifications of the standard least squares procedure has 
been proposed 7 -^. Currently the most common practice 
is to use the Maximum Entropy algorithm for the ana- 
lytic continuation problem 9 . It's based on Bayesian in- 
ference concept and essentially uses a priori knowledge 
about the properties of the spectral function, its positiv- 
ity and sum rules. In other words, the MaxEnt algorithm 
is essentially a nonlinear formulation. Existing realiza- 
tions of MaxEnt do not allow the analytic continuation 
of functions, whose norm is not a known constant (self- 
energy, for instance), or do not have a definite sign (for 
example, off-diagonal components of Green's function). 

The aim of this paper is to introduce a new analytic 
continuation algorithm, which is linear with respect to 
input data and strong enough to be considered as an al- 
ternative for the MaxEnt method. Section II of the paper 
is devoted to the general idea behind the method with- 
out giving details of a particular realization. It explains 
optimal regularization principle, which is built on top 
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of Tikhonov regularization ansatz. Some more specific 
realization details are described in sections III and IV. 
They, in principle, may be revised and adjusted to a cer- 
tain extent for specific problem. In Section V we present 
practical results of analytic continuation for the Hubbard 
model on the Bethe lattice, which is then compared with 
corresponding results from the MaxEnt algorithm. In 
Section VI we conclude the paper. 

II. THE OPTIMAL REGULARIZATION 
MATRIX 

Mathematically speaking, analytic continuation re- 
quires a solution of a linear integral equation for unknown 
function F(u>) 

MF{uj) = F{iuj) (4) 

using certain a priori information. Here M is the ill- 
posed linear operator which analytically continues from 
the real to the imaginary axis, and function F(iu) is ap- 
proximately known from a QMC simulation. 

If we have introduced certain basis of orthogonal func- 
tions to represent F(ui) and F(ioj), the above integral 
equation becomes a set of linear algebraic equations 

Mx = y, (5) 

where M is badly conditioned matrix, and the right hand 
part y is subject to some level of inaccuracy. Formally, 
the rank of the set is infinite. One would expect that for 
a properly chosen basis the effective reduction to a finite 
set is possible. Such a basis, suitable for the analytic- 
continuation problem is introduced in Section III. 

The staring point of the proposed method is a 
Tikhonov regularization of the problem-^. Let us search 
for a vector x, which minimizes the Tikhonov functional: 

T[x;R] = \\Mx - y\\ 2 + (x,Rx) (6) 

Here R is a regularizing Hermitian matrix. Vector y is 
known approximately: y — y + Sy, where y is the mean 
value of the vector y, and the deviation Sy is a random 
quantity distributed with zero mean value and character- 
ized by the covariation matrix K y , such that: 

5y~=0, K y = SyW (?) 

(hereafter the line over an expression denotes the QMC 
expectation) . 

Varying the functional J-[x; R] with respect to x, one 
obtains a condition for the vector x, which gives the min- 
imum of the functional for given y and R: 

x = XM^y, where X = (M^M + R)- 1 (8) 



part: Mx = y. We average the mean square deviation of 
x from x over all possible values of the random vector 5y 
to obtain the following: 

\\x-x\\ 2 = Tr{xAx - 2XB} + Tr{xx^} (9) 
A = M^Mxx^M^M + M^KyM (10) 

B = M^Mxx^ (11) 

It should be obvious, that a proper choice of regular- 
izing ma trix R is r equired to provide a satisfactory small 
value of ||x — x\\ 2 . On the other hand, a desired type 
of the behavior of a solution x, cannot be determined 
or even estimated from the results of QMC simulation 
alone. Construction of a regularizing algorithm presumes 
utilization of certain a priori information about the prop- 
erties of the solution. 

Very generally, the use of a priori information means 
making an assumption that the result x is not arbitrary, 
but falls into certain classes of possible solutions. For 
example, one can suppose that the resulting function is 
smooth, not too large, etc. We denote by (. . .) the aver- 
age over possible class of the solutions. It is reasonable 
to require that the regularizing functional R delivers the 
minimum of the deviation ||x — x\\ 2 in the average over 
the set of possible solutions, 

(||x-5|| 2 ) = Tr(X{A)X - 2X(B}) + Tr{xx^) = mux 

R 

The main idea of the proposed method is to solve this 
minimum problem with respect to R and then to obtain 
x from ||5J). A specific choice of the possible class of solu- 
tions is discussed in Section IV. However already at this 
point it is worth noting that we actually do not need a 
complete information about the possible solutions, but 
only the first and the second momenta of x, since only 
these quantities appear in (Tj"2")) (recall expressions (fTUI 
HP for A and B). 

Solving (|12[) with respect to R and taking into account 
the additional condition 6R = dR? , one obtains an equa- 
tion for the matrix X (the trivial solution X = does 
not relate to the problem): 

(A)X + X{A) = (B) + (B t ) (13) 

In this way we formally found a system of N(N + 
l)/2 linear equations for all elements of the matrix X 
(whereN is the dimension of the square matrix (B)). 
However, there is a more efficient way to solve the ob- 
tained equation (|13() . 

Let us denote eigenvalues of the matrix (A) by Xi and 
perform an unitary transformation to a primed basis, in 
which (A) is diagonal. In this basis a solution of equation 
(fT3| can be expressed explicitly: 



Assume that the vector x is an exact solution of the 
set of equations |5]| with an exactly known right-hand 
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Changing to the eigenbasis of the matrix (A) allows us 
to organize the numerical solving of the system (|13|) in an 
efficient way. The search procedure for eigenvalues and 
eigenvectors of (A) is quite stable because (A) is Hermi- 
tian. The advantage of such an approach in comparison 
with the direct solving of (fT3|) is that instead of solving 
a system of N(N + l)/2 equations (oc N 6 operations) it 
is enough to diagonalize the N x N matrix (A) (oc TV 3 
operations). 

We should note that, mathematically, formula (fT4|) un- 
ambiguously gives finite expressions for all the elements 
of X, since its denominator is positive. Indeed the matrix 
M'M is nonnegative definite. Matrix (xx') is a corre- 
lation matrix, which is also essentially positive definite. 
Finally, the matrix K y is also positive definite, if the vec- 
tor y is defined with nonzero inaccuracy (all components 
having nonzero dispersion). In this way we ensure that 
the matrix (A) is positively defined and therefore all its 
eigenvalues are positive. 

On the other hand, the practical realization of the 
method faces certain difficulties, because the numerical 
solution suffers from round-off errors in our calculations. 
The singular value decomposition procedure does not 
help much. To get rid of these round-off errors, we have 
switched to performing the calculations with to arbitrary 
number of digits, using CLN library version 1.1.131*. 



III. CHOICE OF REPRESENTATION. 
CORRELATION MATRIX 

Now we apply results of the previous Section directly 
to the analytic continuation of a function F(u>) from the 
imaginary axis to the real one, assuming F(u>) to be an- 
alytic in the upper plane. QMC simulation gives values 
of the function F at Matsubara frequencies iu>k on the 
imaginary axis: 



F k = F(iuj k ), 



(15) 



Let Wo be a typical frequency scale of the problem. Let 
us introduce a conformal mapping of the upper frequency 
plane to a circle of a unitary radius with the center at 
zero point: 

UJ-IUJQ . l + z 

oj — > z : z = . ijj = iojq- (lb) 

CO + IUIq 1 — Z 

The mapping is illustrated by figure [T] In this case 
all imaginary frequencies are mapped to a segment z 6 
[— 1 ; 1] ; all real frequencies correspond to the circle with 
radius 1. 

If F is analytic in the upper plane as a function of 
complex frequency, it has the same property in the unit 
circle in the z-plane and can be expanded in the Taylor 
series at point z = 0: 



F(uj( Z )) = J2fnZ n , fn = 
n=0 



1 

27ri 



F(u(z)) 



dz (17) 




FIG. 1: Schematic view of the conformal mapping of upper 
frequency plane to a circle. 



Assuming that the function F(z) is smooth, we can 
then take into account a finite number of terms in this 
expansion. The proper number TV of terms to be kept 
can be estimated from test runs of the program for a 
known function F. After the expansion coefficients f n 
are determined, we can then sum up the series at any 
point of the circle \z\ = 1 and thus restore values of F(ui) 
on the real axis. 

The chosen representation M entering equation (JSJ) is 
then a K x N matrix with elements given by formula 



Mkn = 



( UJ k - U)Q 

\u!k + ujo 



(18) 



Using some a priori information about the expected 
solutions F(u> £ K), we can then choose a correlation 
matrix (F(uj)F*(u;')) | w , W 'eR - This correlation matrix 
also can be expanded into the double Taylor series at 
zero point, similarly to (jTTJ) : 



(F(u,)F>')> 



oj,oj' 'GIF 



N-l 

E 



{f n f*,)z n ( z 'y 



(19) 



1 



(2ttz) 2 



(F(w(z)) J F 1 *(w(z')))z" (n+1) z' ( "'^ 1) dzdz' 



|*|,|*'|=i 



(20) 



\z\ = l 



Practically, it was found that w 30 — 50 terms is 
enough for the calculation. A value of ojq was not found 
to be really important; it is sufficient just to take loq of 
the right order of magnitude. 
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IV. CORRELATION MATRIX OF 
LORENTZIAN PEAKS 

The explicit form of the correlator (F(uj)F*(lu')), lead- 
ing to satisfactory results, can be obtained from the fol- 
lowing procedure. 

Let us assume that function F is a superposition of 
several Lorcntzian peaks having the same width 7 
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The positions of the peaks Clj and their magnitudes Z\ 
are random quantities with known statistical properties 



<F( Wl )F*(w 2 )> 



j 

E 



jjti n i +il^ 2 - Qj' - il 1 ZSI 

(22) 

We assume that all Zj are equally distributed and un- 
correlated with each other, so that two model parame- 
ters (Z 2 ) and (Z) completely determine the distribution 



(21) of Z 3 



^— ' \ oj\ — Qj + i-y 0J2 — f2i — il 
j=i 



<^> 2 E 



1 



OJl — Qj + J7 Ll>2 — 



vy 



(23) 



Finally we assume, that the values flj are distributed independently with certain model distribution, for example 
the Lorentzian: 



1 



1 



nQ M 1 + (nj/n M y 



(24) 



where Om is the scale of the spectrum. For example, half of the Hubbard U is a good guess for this quantity in 
Hubbard-like models. 

Combination of formulas (p?0"|) and (|23[) yields a result suitable for a practical calculations: 



(fnf*,) = J(Z 2 ) 



P(n)/(n,7,Q)rV,7,ft) dn+ 

/+00 r + cc 

P(ni)/(n,7,ni) dOi / P(fi 2 )I«7,Q 2 ) dfl 2 (25) 
-OO «/ —OO 



(»7-q- 



-2iojQ 



71 = 

n ^ 



Although these integrals could be calculated analyti- 
cally for a particular Lorentzian form of P(£l), the an- 
swer is very complicated and has significant computa- 
tional complexity. It is therefore much more practical to 
perform an integration numerically for every pair of n 
and n' . 



was used as an impurity-problem solver. This version 
of the QMC solver is particularly suitable, because it 
can produce high-accuracy data for Q(iuj n ) directly in 
iw-domain, up to quite high Matsubara frequencies. 

The regularization algorithm uses the model described 
in the previous section with the Gaussian distribution of 
flj . An important peculiarity of the approach is that we 
work not with the Green's function, but with the quantity 
E — A, defined by the relation 



V. MOTT TRANSITION AT BETHE LATTICE: 
PRACTICAL CALCULATION OF DOS 

To illustrate the potential of the method, we recon- 
struct the density of states for the half-filled single-band 
Hubbard model defined on the Bethe lattice. Near- 
neighbor hopping constant for the model was equal to 
0.5, and the Hubbard U ranged from 1.0 to 3.0. Inverse 
temperature in the simulations was j3 = 50. The Mott 
transition was observed between U = 2.4 and U = 2.5. 
Data for Q(t) and Q(iui n ) was produced by the DMFT 
self-consistent loop. The continuous-time QMC code 12 



iuj n + A(iu n ) - E(iw n ) 



(26) 



Practically such a method has an advantage over the di- 
rect continuation of Q(iu} n ), because E — A does not have 
trivial l/(iw) asymptotics. Therefore, a smaller number 
of z-harmonics is required for an adequate representa- 
tion of this function. This also provides us with a phys- 
ically important property J^ oQ A(uj)duj — 1, assuming 
that A — E remains finite as lu — > 00. Figure [2] illustrates 
the method, and Figure [3] shows the results for several 
values of U corresponding to metallic behavior. Unfor- 
tunately, the procedure works well only for the metal 
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(b) 

FIG. 2: (a) Input data for T,(iuj n ) - A(ico n ) at U = 2.3. (b) 
Reconstructed real and imaginary parts of E(o») — A(w). 300 
Matsubara frequencies were used, with an estimated error- 
bar in E — A is about 10 -2 , the noise at different Matsubara 
frequencies was assumed to be uncorrelated. In the regular- 
ization procedure, the eigenvalues of 40 x 40 matrix have been 
obtained, using an accuracy to 80 decimal digits. 



phase, because the insulator gap in DOS requires an in- 
finite S — A, which is clearly an impossible task for the 
continuation procedure. 

For comparison, we present the result of MaxEnt cal- 
culation for the same QMC data. We have used a rel- 
atively simple MaxEnt program^. This code does not 
diagonalize the covariance matrix and uses a full search 
algorithm. 

As can be seen from Figure 3, our procedure produces 
the data with much more details. Firstly, in all graphs 
the DOS decreases more rapidly at large energies. This 
is quite physical, one would expect that energy bands 
should be quite well confined. Secondly, in our case 
the shape of the dips at small energy in Figure [3] b,c is 
more similar to a "pre- formed" insulator gap. Finally, 
we were able to resolve very pronounced DOS shoul- 
ders near the onsets of the Hubbard bands. This agrees 
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FIG. 3: Reconstructed densities of states for U = 1.0, U = 2.2 
and U = 2.3. Parameters of the regularization procedure 
are the same as in Fig. 2. Results of the regularization are 
represented by solid lines. MaxEnt predictions are shown with 
dots. 



well the with the predictions of dynamical density-matrix 
renormalization-group theory**. 

Of course, it should be noted that the particular ap- 
pearance of the resulting graphs depends on a choice of 
a priori parameters Qm & n d 7- However, all qualita- 
tive peculiarities discussed in the previous paragraph are 
found to be quite robust against the variation of & n d 
7 in a reasonable interval (approximately, from 0.1 to 
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2). Values of ujq, (Z) 2 and (Z 2 ) do not affect the result 
considerably, and so we put these quantities to 1 in the 
calculations. 

VI. DISCUSSION AND CONCLUSIONS 

To conclude, we developed a new regularization ap- 
proach for the analytic continuation of numerical data 
from imaginary to real axis. The method is the best pos- 
sible regularization in sense of (| 1 2[) . so that we produce 
an optimal regularization matrix given a class of possible 
solutions and errorbar of the input data. 

An important property of the optimal regularization 
approach is its linearity with respect to the input data. 
This makes possible the continuation of the self-energy 



and off-diagonal Green's function components. In princi- 
ple, we could follow MaxEnt paradigm and put additional 
constrains like the requirement A{u>) > 0. However this 
would make our scheme nonlinear and seriously reduce 
its flexibility. There is also negative consequences of the 
linearity of our approach: unphysical regions of negative 
DOS can occur for certain input data. It should also be 
noted that, accordingly to our observations, our method 
requires higher accuracy and larger number of Matsubara 
frequencies. In this case, our scheme has an advantage, as 
it is illustrated by the practical calculation for the Mott 
transition. 
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